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This paper describes a mixed direct-iterative method for boundary integral for- 
mulations of dielectric solvation models. We give an example for which a direct 
solution at thermal accuracy is nontrivial and for which Gauss-Seidel iteration 
diverges in rare but reproducible cases. This difficulty is analyzed by obtaining 
the eigenvalues and the spectral radius of the iteration matrix. This establishes 
that the nonconvergence is due to inaccuracies of the asymptotic approximations 
for the matrix elements for accidentally close boundary element pairs on different 
spheres. This difficulty is cured by checking for boundary element pairs closer than 
the typical spatial extent of the boundary elements and for those pairs perform- 
ing an 'in-line' Monte Carlo integration to evaluate the required matrix elements. 
This difficulty are not expected and have not been observed when only a direct 
solution is sought. Finally, we give an example application of these methods to 
deprotonation of monosilicic acid in water. 



1 Introduction 

An interesting development in computational molecular biophysics over the 
past decade has been the surprising utility of dielectric models of solvation of 
molecular solutes in waterErt3 This is surprising a priori because this approach 
neglects almost all of the molecular characteristics of solvation. A posteriori 
molecular calculations hava^b£come available checking the basic soundness of 
the dielectrkimDdel result£j"0 and checking features of the underlying molec- 
ular theoryExEj 

Arguments that support such models are simple and broad: much of the 
solvation phenomena in water are dominated by electrostatic interactions. 
These models provide a physical description of solvation of electrostatic in- 
teractions. If we permit a macroscopic empirical parameterization then they 
are indeed useful. Furthermore, dielectric models permit a conceptually nat- 
ural and feasible coupling of solvation theory with electronic structure tools 
of traditional computational chemistry. For these reasons too the dielectric 
models have been helpfulcHi 

The numerical challenge in applying these models is the solution of the 
Poisson equation 

V»e(r)V$(r) = - 47rp(r) (1) 

where p(r) is the density of electric charge associated with the solute molecule, 
e(r) gives the local value of the dielectric constant, and $(r) is the electric 
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potential. This equation can be a challenge because the function £(r) changes 
abruptly on the modeled molecular surface of the solute and that surface must 
sometimes exhibit nontrivial variation on an atomic scale. 

Because the most important difficulty is associated with treatnifint of the 
molecular surface, boundary integral methods are advantageousoE3 Those 
methods permit the concentration of numerical resources on the description 
of the molecular surface. The resolution in the description of the molecular 
surface can then be directly associated with the accuracy of the numerical 
calculation. 

The accuracy requirements of relevance to us are associated with confor- 
mation free energy differences comparable to k^T and with treatment of the 
effects of molecular solvent structure hjt. integrating out probe water molecules 
with the help of this dielectric modeleil These interests put high demands of 
accuracy and speed on the numerical methods. Stringent testing of the accu- 
racy of these dielectric_models for structural optimization has been pursued 
only relatively recentlycil 

It might be questioned whether it makes sense to solve approximate di- 
electric models to the accuracy discussed here. We offer two responses. First, 
though the model is approximate, attempts to draw conclusions from the model 
results are complicated by non-physical errors superposed on the model results. 
Second, if the model results are valid enough to be helpful, then they might 
serve as^n initial approximation upon which more refined treatments might 
be built Cj In that case, understanding the accuracy of the initial predictions 
would be important. 

The accuracy that can be achieved in the solution of Eq. (|l|) through a 
direct boundary integral approach will be limited by the dimension of the set 
of linear equations that corresponds to the linear boundary integral equation. 
Since the dimension of that set will be relatively small if direct solution methods 
are used, substantial numerical resources can be invested in obtaining accurate 
coefficients for the linear equations. 

When direct methods become unfeasible, it is natural to apply iterative 
solution methods to the direct solution used as an initial estimate. Because 
the initial solution is expected to be good, the iterative effort is expected to 
be modest. Iterative methods permit a larger number of linear equations. But 
numerical sophistication in the evaluation of the larger number of matrix ele- 
ments becomes prohibitive. Thus, the price to be paid for the higher resolution 
is that the most matrix elements are obtained 'on the fly.' The methodological 
problem of this paper is the formulation of the mixed direct-iterative methods 
for solution of Eq. (|l]) ; and the identification and correction of a difficulty that 
can arise. 
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Figure 1: Ca^^ ■ ■ -CI" potential of mean force in water at normal pressure and 298K. The 
curve labeled 'molecular dynamics' is the literature molecular resulted! The curve labeled 
'dielectric' is obtained by the revised method of Section 2.4. 



The results of Fig. |l| present an example that will be used because an it- 
erative difficulty can be reproducibly exhibited. Shown there are calculations 
of the potential of mean force between a Ca"''^pipn and a Cl~ ion in waterEa 
These results utilize the van der Waals surfacecll and the radii recommended 
by Rashin and HonigE3 We include here some qualitative notes about physical 
aspects of these results. Firstly, we have much less experience with simulation 
results for this potential of mean force than we do, for example, with Na"*' • • • 
Cl^. Thus, we view the simulation results of Fig. as preliminary. Assuming 
those results are born-out by further study, the interpretation would be that 
the Cad:+ holds its solvation shell sufficiently tightly that no contact minimum 
existsLJ Secondly, the dielectric model result predicts an over deep contact 
minimiirn Tti tViis rcspcct the present results are consistent with previous com- 
parisonsOEjOo and these results are therefore not newly troubling. Thirdly, 
the maximum in the dielectric model result near 3.8A where the spheres just 
touch is expected to be correct though it is clearly a subtle feature on the global 
scale shown here. The free energy of the separated ion pair is approximately 
700 ksT so resolution of such features here requires an relative accuracy of 
about 0.1%. Still, the relative height of that maximum is not negligible on a 
k^T energy scale. Calculations that establish the correctness of k^T features 
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require careoSlflJ Such methods are the topic of this paper. 
2 Methods 

Here we catalog the methodological results used below for solution of the Pois- 
son Eq. (01. Further discussion of the genesis of these results can be found 
elsewhereEll We first cast Eq. ™) as an integral equation, e.g.: 



e(r)$(r) = $(")(r) + j 



(r-r')»V'e(r') 



47r |r - 



$(r')dV'. (2) 



The quantity <i>*^°^ (r) is the electrostatic potential in the absense of the medium. 
Because the model assumes that e(r) has a sharp step at the molecular surface 
the integration on the right collapses to a 2-dimensional integration over the 
molecular surface. That molecular surface is defined as the boundary between 
the molecular volume — modeled as the union of spherical volumes centered 
on solute atoms — and the solution region. For r infinitismally outside the 
molecular surface, Eq. (||) provides a closed equation for <l>(r) on the molecular 
surface. Once $(r) is obtained on the molecular surface, it can be used on the 
right side of Eq. (g) to construct the potential elsewhere. 

From such solutions we construct the interaction part of the chemical po- 
tential of the solute as 

A/i^") = (0 / P(r) ~ <i>„(r)) dV. (3) 

The subscripts I and v indicate 'liquid' and 'vapor,' respectively, so that this 
difference is the electric work required to charge the solute in the liquid relative 
to the vapor. This requires the solution of Eq. (j^) twice, once for the liquid 
with 

ei{T) = em + {ss - em)ri{y), (4) 
and once for the vapor with 

e«(r) = Em + (1 - em)?7(r). (5) 

Here ry(r) is a step function that is one outside the molecular volume and zero 
otherwise; Eg is the dielectric constant of the solution and e™ is an assigned 
'dielectric constant of the molecule.' The latter parameter is used to match 
the polarizability of the solute. The formulation-Eq. (^) makes it simple to 
match a given polarizability by adjustment of This is because the kernel 
is proportional to the electrostatic potential due to a surface dipole density. 
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Thus, we perform calculations with £u(r) and with $'^*'^(r) chosen to describe 
a uniform external electric field. The induced electrostatic potential in the 
far field is associated with the induced dipole moment. The correlation of the 
induced dipole moment with the external field strength provides the modeled 
molecular polarizability. 

2. 1 Rules for coarse direct calculations 
A discretized version of Eq. (||) is 

e,$(r„) = $(")(r,)+^«;«;3$(r;3). (6) 

Here Tq, is the a-th 'plaque point' — a point on the molecular surface obtained 
by a uniform sampling, for exampjjahi exploiting either quasi-random number 
series or 'good lattice' procedures The plaques are defined as the Voronoi 
polyhedra of the plaque points on the nonburied surface of each sphere. The 
matrix of coefficients Wq,^ can be obtained as follows: 



R{sp)'^{es - Em) (r^ - r^) • 

M{sp) 



«^"/3 = TTTT-^ 1^ ^^.13 ' "^Z^' (^) 



{ie/3} 

and 



f^TT^f) Ed— s^»)-^/^ (8) 



2V2 V M{sa) 



These are Monte Carlo estimates of plaque integrations. Further details can 
be found elsewhereEil The set of sampling points that are within the plaque /? 
is denoted by {i S /3}. M{sg) is the number of points on the sphere s^g that 
supports plaque /3. In Eq. (j^), dia is the angle between the surface normals 
at plaque point a and the sampled point i. That formula arranges to use 
sample points outside the plaque to calculate the solid angle subtended by 
that plaque at the plaque point. The purpose is to reduce the variance of the 
Monte Carlo estimate. All sample points on the surface of the sphere that 
supports plaque a should be used in the estimation. But whether any sample 
point resides on plaque a depends on resolution of buried surface because the 
plaque boundaries sometimes follow the boundaries between the exposed and 
buried surface. 

For accurate calculations on small molecule solutes, we have found that 
the computational time is dominated by the Monte Carlo effort. Thus, we 
use these formulae differently to avoid some of that effort. We use Eq. (||) 
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with no Monte Carlo sample in addition to the plaque point; these formulae 
constitute one point estimates then. However, the set of plaque points is 
now a larger set of equidistributed surface points, elements of a fine lattice. 
We form the equations to be analyzed by contraction of that description to 
one based upon coarse plaques constructed from an initial sequence of the 
fine lattice points. We then require the equality of the potential at all fine 
lattice points residing on the same coarse plaque. This requirement results 
in an overdetermined system. This system is analyzed with a singular value 
decomposition to obtain tha^plaque potentials that minimize the mean square 
residual of those equationsca 

The advantages of this approach are that much of the Monte Carlo ef- 
fort can be avoided and that some account is taken of the spatial variation 
of $'^''^(r) within a coarse plaque. The principal disadvantage is that this 
approach requires more memory and this disadvantage can be severe. 



2.2 Rules for Gauss-Seidel iteration 

Here we give the rules used in our iterative calculations. We begin with an 
approximate solution $(ra). That approximation is then updated in place 
according to 

$(r„) ^ (es - w^^y' J $(")(r„) + Wc^'i'irp) I (9) 

sequentially for all a. In this calculation the plaque points are the points of 
the fine lattice. Because that set of points is expected to be large, the off- 
diagonal coefficients Waf3 are evaluated 'on the fly' using Eq. ^ but no Monte 
Carlo sample in addition to the plaque point. Our experience is that accurate 
evaluation of the diagonal coefficients Waa is important, so we evaluate them 
at an initial stage of the calculation and store them for later use. 

These methods were applied to the calculation of the pair potential of 
the mean forces between a Ca++ and a Cl^ in water. It was found that 
the Gauss-Seidel iteration scheme converged almost always, but diverged in 
rare but reproducible circumstances. A necessary and sufficient condition for 
the convergence of Gauss-Seidel iteratiDn is that the spectral radius of the 
iteration matrix must be less than onecj Plotted in Fig. || are the eigenvalues 
of the Gauss-Seidel iteration matrix for the Ca"'""'" • • • Cl~ problem for the case 
of r = 3.7A with 36 plaque points. One eigenvalue is much greater than one. 
Also plotted there are eigenvalues of the Gauss-Seidel iteration matrix for the 
same circumstances except that matrix elements were obtained by a modified 
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method described below. All eigenvalues are now substantially less than one. 
Using the corrected methods the divergences have not been observed. 



2.3 What we did to fix the problem 

It was when those Monte Carlo efforts were economized that iterative diver- 
gence occasionally presented a problem. Furthermore, the diagonal matrix 
elements are always calculated the same way. This suggests that the observed 
difficulty was due to the one-point estimate of the off-diagonal elements used 
when the iterative calculation was implemented. 

The estimate Eq. will have a larger variance the closer the point to 
plaque (i. Thus, it is reasonable to suspect those matrix elements corresponding 
to close aj3. It was verified that this suspicion is correct by replacing Wafi 
by "Wfjp whenever Tq, is on a different sphere than plaque (3 and \Ta — r^| < 
2i?(s^)/ ■y/Af(s^). This is a statistical estimate of the radial extent of plaques 
on center S/3. This unsatisfying maneuver eliminated the divergence. 

A geniune solution is to implement a Monte Carlo calculation of those 
matrix elements Waf3 identified as potentially problematic. For Tq, close to 
plaque /3, an approach like that of Eq. using points sampled outside the 
plaque would be most appropriate. But for far from plaque /3, it would be 
natural to use points on the plaque. That we could use either approach when 
the point Va is not on the plaque is justified by the relation 



(r-rQ. V'£(rO 
47r|r- r'l^ 



d^r' = 0, 



(10) 



valid for r outside the sphere. This is an application of Gauss's law. Thus we 
could estimate the required integral using either points on the plaque (} or on 
a complementary spherical surface. Ilsing 1/2 of each estimate is an example 
of the method of antithetic variates:! 



a/3 



R{spf{es - £m) 
2M{sp) 



E 



r„ - r,: 



E 



(11) 



_{ie/3} ' " " {tm ' " 

This effort is expended only when is on a different sphere than plaque /3 
and Ire, - r^| < 2R{sp)/ y^M{Ipj. 
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Figure 2: Upper panel: Eigenvalues of the Gauss-Siedel iteration matrix obtained as de- 
scribed in Section 2.2 for 36 plaques on the Ca^^ ■ ■ • Cl~ di-ion for r = 3.7A. One eigen- 
value is much further from the origin than 1.0. Lower panel: Eigenvalues of the Gauss-Siedel 
iteration matrix obtained by the revised approach described in Section 2.3. 
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3 Deprotonation of monosilicic acid in water 



As an application of the methods developed in the previous sections we will 
treat the deprotonation of monosilicic acid in water at 298K 

Si{OH)i^ Si{0H)-30- +H+. (12) 

The monosilicic acid molecule and anion are depicted in Fig. ^ The equilibrium 
ratio is 

[Sz{0H)30-] [H+] 

with concentrations in molar units. The quantity we seek is the free energy of 
reaction 

AG(°) = -RTlnK (14) 

measured to be 13.5 kcal/mol0'@ 

This is a helpful example for several reasons. Acid-base equilibria are 
an important application of these models in molecular biophysics. Addition- 
ally, these solutes have not be treated previously by these methods. Thus, 
the expectations for the radii-parameters required can be tested outside the 
conventional parameterization suite of solutes. 

3. 1 Solution thermodynamic formulation 

It is more physicalEl to consider the reaction 

Si{OH)i + H2O ^ Si{0H)30- + H3O+. (15) 

The reaction described this way does not result in a net loss of chemical bonds 
and this is likely the case in waterEil The ratio 

- ^ [SijOHhO-] [H3O+] 

[St{OH),] [H2O] ^ °^ 

is then dimensionless. This equilibrium ratio may be obtained afl 

k{T) ^ i^(°)(T)exp [-AA/^(^V-Rrj , (17) 

where K^'^'' (T) is the ideal gas result obtainable from standard formulae^ and 
AA^(-) ^ A^g^, + A^,ill^^^^_ - A^^^Sl^ - A^^i%^^^ . (18) 

K and K are related by if = K/[H20] and the free energy of reaction is simply 

AG(°) = -RTlnK - RT\n[H20]. (19) 

We assume that the solute concentrations are sufficiently low that the formal 
concentration of H2O is satisfactory. 



Table 1: Partial charges for Si(OH)4 and Si{0H)30-. 



Atom 


Si(0H)4 


Si(0H)30- 


Si 


1.62 


1.52 


01 


-0.90 


-0.89 


02 


-0.90 


-0.90 


03 


-0.90 


-0.92 


04 


-0.90 


-1.06 


HI 


0.49 


0.44 


H2 


0.49 


0.41 


H3 


0.49 


0.41 


H4 


0.50 





3.2 Electronic structure results on the isolated molecules 

All electronic, structure calculations were performed using the GAUSSIAN- 
92 programE3 Two different quantum mechanical methods were employed: 
Hartree-Fock (HF), and HF followed by a second-order order MoUer-Plesset 
(MP2) correlation energy correction. Two different basis sets were used: 6- 
31G(d) [also denoted 6-31G*], and 6-31G+-H(2d). The "(d)" and "(2d)" de- 
notes tliat the 6-31G basis is supplemented by one and two sets of polar- 
izationE3 d-functions, respectively, on the heavy (non-hsdrogen) atoms. The 
"++" denotes that the basis is supplemented by diffusecJ functions. Together 
the quantuin mechanical method and basis set specifies a theoretical model, 
e.g., SauerEa has performed HF/6-31G(d) calculations on monosilicic acid. The 
optimized geometries for H2O, H3O+, Si(0H)4 and Si(0H)30^ were deter- 
mined by analytic gradient techniques using the HF/6-31G(d) model. The 
structures found for monosilicic acid and its anion are shown in Fig. |^. The 
bond distance and angle faLH20 is 0.947A and 105.50°, respectively compared 
to the experimental valuesEj of 0.957A and 104.5°, respectively. The calculated 
value for all three H-O-H bond angles for H3O+ is 113.06°. Teppen et aW3 
have examined the effects of basis set size and electron correlatiDn corrections 
on the properties of monosilicic acid. For a 6-31G(d) basis, theyl£3 find that the 
MP2 bond lengths for Si-O and O-H were 0.022 and 0.023A larger, respectively 
than the HF values and the ME2 Si-O-H angle decreased 2.9° from ±he HF 
value. For the HF method, theyS also find that the MG6-311G(2d,2p)0 bond 
lengths for Si-O and O-H were 0.007 and 0.009A smaller than the 6-31G(d) 
values and the MG6-311G(2d,2p) Si-O-H bond angle increased 1.3° from the 
6-31G(d) value. For the present work these differences are acceptable, and 
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Table 2: Partial charges for H2O and H3O+. 



Atom 


H2O 


H3O+ 





-0.82 


-0.62 


H 


0.41 


0.54 



Table 3: Electronic energies and polarizabilities. 



Molecule 


Eq (hartree) 




Si(0H)4 


-590.89169 


6.9 


Si(0H)30- 


-590.29838 


7.4 


H2O 


-76.01075 


1.44 


H3O+ 


-76.28934 


1.44 



HF/6-31G(d) was used to optimize geometries. 

To obtain atom-centered charges (Tables ^ and ^ necessary for the sol- 
vation enexgy calculation, a fit of the HF/6-31G(d) electrostatic potential 
(CHELPGE3) was performed. Harmonic vibrational frequencies and rota- 
tional constants were computed with the HF/6-31G(d) model and were used 
to compute the partition functions in Eq. (p7|). Since the HE na£.thod over- 
estimates frequencies, the computed values were scaled by 0.88B The elec- 
tronic ground state energies, Eg, (Table ^ were calculated by optimizing the 
molecules with the MP2/6-31G(d) model. The calculation of polarizability, 
a = {oixx + Q!yy -t- azz)/3, is sensitive to the basis set. Eor H2O at the HF/6- 
31G(d) geometry, a = 0.70, 0.90, 0.90, and, I.OGA^, calculated with the 6- 
31G(d), 6-31-hG(d), 6-31G(2d), and 6-31-h-hG(2d) basis sets, respectively The 
6-31-|--|-G(2d) value agrees reasoiiahli' well with another HE calculationfj a = 
1.17A3. The experimental valueH'El for H2O, a = 1.44A3, was used in the 
solvation calculations for both H2O and H3O+. Eor Si(0H)4 and Si(0H)30~, 
a was calculated with the HE/6-31G+-t-(2d) model using the HE/6-31G(d) 
geometries. These values were scaled by 1.36, the ratio of the experimental 
and HE/6-31G-h-h(2d) a values for H2O. The scaled values (Table |) were then 
used in the solvation calculations. 



3.3 Results for the deprotonation of monosilicic acid 

Three different calculations were performed for the free energy of deprotona- 
tion, Eq. (12). The first calculation did not include molecular polarizability. 
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Figure 3: Monosilicic acid molecule (upper) and anion (lower) established by the electronic 
structure calculations of Section 3.2 
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nor spheres on the H atoms of Si(0H)4 and Si(0H)30^. The second calcula- 
tion included molecular polarizability, but again not spheres on the H atoms. 
The final calculation included both molecular polarizability and spheres on 
every atom of Si(0H)4 and Si(0H)30~. This sequence of calculations reflects 
our chronological approach to this system starting with the simplest model, 
gradually adding more complications, and gradually refining the values of the 
parameters used. This approach also gives helpful information on the sensitiv- 
ity of the calculation to the empirical parameters used. 

The evaluation of the vibrational, rotational, and translational partition 
functions of the isolated molecules on the basis of the electronic structure 
results leads to a multiplicative contribution of 1.30 to the equilibrium ra- 
tio of Eq. ([TtI). The change in the electronic energy AEq was found to be 
197.5 kcal/mol. 

In all of the calculations the water and hydronium ions were treated as 
single spheres of radius 1.6 A on the O atom. A molecular dielectric omatpit Sm 
was 2.42. This reproduces the polarizability of the water moleculeoEJ'E£l The 
solution dielectric constant was Es = 77 A appropriate to water at 298K and 
l.Og/cm'^. The calculation of the excess chemical potential of solvent species 
used 186 coarse lattice points and 936 fine lattice points. Ten Gauss-Seidel 
iteration passes were applied to the coarse solution in all calculations. The 
difference in excess chemical potential between the hydronium ion and water 
was found to be -107 kcal/mol. All calculations on Si(0H)4 and Si(0H)30^ 
used 36 coarse lattice points and 936 fine lattice points on each sphere. 

The first calculations for Si(0H)4 and Si(0H)30~ used a sphere of radius 
l.SA on the each Si atom and a sphere of radius 1.65A on each O atom. 
Em =1.0 was adopted, thus ignoring the polarizability of the molecule. The 
difference in excess chemical potential was found to be -49.6 kcal/mol. This 
calculation gives 38.3 kcal/mol for the change in free energy for Eq. (|lj), in 
poor agreement with experiment. 

In the next calculation a sphere of radius 1.8 A was centered on the each 
Si atom and a sphere of radius I.60A on each O atom. Em =2.95 and 3.40 
were assigned to the Si(0H)4 and Si(0H)30~ molecules, respectively. These 
values match the estimated molecular polarizabilities discussed in the previous 
section. The difference in excess chemical potential was found to be -57.9 
kcal/mol, which gives 30.0 kcal/mol for the change in free energy for Eq. (p^. 
This improves upon the calculation which ignored the polarizability of Si(0II)4 
and Si(0II)30^ but still differs from experiment by more than a factor of two. 

The final calculation had spheres on every atom of the solute molecules. 
The O atoms were given radii of l.lA, the Si atoms were given radii of 1.8A, 
and the H atoms were given radii of 1.3A. Values of 3.10 and 3.55 were used for 
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the Sm of Si(0H)4 and Si(0H)30^, respectively. Again these values were de- 
termined to match the polarizability of the molecules. The difference in excess 
chemical potential between Si(0H)4 and Si(0H)30^ was -68.1 kcal/mol. This 
calculation gives 19.8 kcal/mol for the change in free energy for the deproto- 
nation of Si(0H)4 in water, a much improved agreement with experiment. 

4 Conclusions 

The iterative divergence occasionally encountered in the calculation of the 
Ca+^ • • • Cl^ pair potential of mean force was due to inaccuracies of the asymp- 
totic approximations used for the matrix elements for accidentally close bound- 
ary elements on different atomic spheres. This problem is cured by checking for 
boundary element pairs closer than the typical spatial extent of the boundary 
elements and for those pairs performing an 'in-line' Monte Carlo integration 
to evaluate the required matrix elements. These difficulties are not expected 
and have not been observed when only a direct solution is considered. 

These methods can give a reasonable description of the free energetics of 
the deprotonation of monosilicic acid in water. A modeled solute polarizability 
and spheres on hydroxyl protons have been found to be important in achieving 
a reasonable agreement between model and experiment. The discrepancy re- 
maining provides a suggestipn of more specific solute-solvent interaactions. It 
has been noted previouslyC^I how specific molecular solvation structure can be 
reintroduced into these models. Those ideas for integrating-out solvent degrees 
of freedom suggest the electronic structure calculations should be preformed 
on complexes of the solutes of interest plus a probe water moleculecil Those 
approaches will require a substantially larger computational effort. 

The necessity of better treatment of the molecular solvation structure is 
also clear in the example of Ca++ • • • Cl^. The probe water molecule approach 
mentioned above would help here too but would require accurate, rapid calcu- 
lations on larger, more complicated solution complexes. It is hoped that the 
methods developed here will make such calculations feasible. 
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